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ABSTRACT 



Magnetic field in the solar corona is usually extrapolated from pliotosplieric 
vector magnetogram using a nonlinear force-free field (NLFFF) model. NLFFF 
extrapolation needs a considerable effort to be devoted for its numerical realiza- 
tion. In this paper we present a new implementation of the magnetohydrodynam- 
ics (MHD)-relaxation method for NLFFF extrapolation. The magneto-frictional 
approach which is introduced for speeding the relaxation of the MHD system 
is novelly realized by the spacetime conservation-element and solution-element 
(CESE) scheme. A magnetic field splitting method is used to further improve 
the computational accuracy. The bottom boundary condition is prescribed by 
changing the transverse field incrementally to match the magnetogram, and all 
other artificial boundaries of the computational box are simply fixed. We exam- 
ine the code by two types of NLFFF benchmark tests, the Low & Lou (1990) 
semi-analytic force-free solutions and a more realistic solar-like case constructed 
by van Ballegooijen et al. (2007). The results show that our implementation are 
successful and versatile for extrapolations of either the relatively simple cases 
or the rather complex cases which need significant rebuilding of the magnetic 
topology, e.g., a flux rope. We also compute a suite of metrics to quantitatively 
analyze the results and demonstrate that the performance of our code in extrap- 
olation accuracy basically reaches the same level of the present best-performing 
code, e.g., that developed by Wiegelmann (2004). 



Subject headings: Magnetic fields; Magnetohydrodynamics (MHD); Methods: 
numerical; Sun: corona 
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Introduction 



The magnetic field configuration is essential for us to understand the solar explosive 
phenomena, such as flares and coronal mass ejections. Besides, the magnetic field also 
plays a crucial role in determining the slowly-evolving structures of the solar corona, such 
as the coronal streamers and the coronal holes. However, direct measurements of these 
magnetic fields are very difficult to implement, and the present observations for the magnetic 
fields based on the spectropolarimetric method (the Zeeman effect and the Hanle effect) 
are basically restricted on the visible surface layer, i.e., the photosphere. Even a routine 
recording of the full surface fields on the photosphere are only available for the line-of-sight 
(LoS) component (e.g., the daily disk magnetogram provided by SOHO/MDI). Most of the 
vector magnetograms at present are recorded locally for active regions and some of them 
may be unreliable because of large random errors and the 180° ambiguity. In view of these 
limitations, researchers hence resort to using physical models to extrapolate (or reconstruct) 
the coronal fields from the observable photospheric magnetogram ( |Sakurai|[l989| [Al^p89 



Amari et al. 1997 McClymont et al. 1997 Aschwanden 2005; Wiegelmann 2008). 



On large scale with relatively low resolutions, the corona fields are usually extrapolated 
from the LoS magnetogram with models including the potential field source surface model 
dAltschuler fc Newkirk|l"969||Hoeksema|1984[ ) and the MHD models ( |Mikic et al.|1999[ [Linker 
et al.|199"9 Feng et al.|2007 2010). By these models and the global map of photospheric field, 
i.e., the synoptic map as usually called, the extrapolated global fields can be used to study 
the general structures of the corona and the heliosphere (e.g., the locations, shapes and sizes 
of coronal holes, coronal streamers and heliospheric current sheet and their evolutions). On 
the local scale with high resolutions, when one's interest is focused on the active regions, a 
most common and powerful way of reconstructing the magnetic fields is the nonlinear force- 
free field (NLFFF) extrapolation from the vector magnetogram. The force-free assumption is 
a good approximation for fields in the low corona but above the photosphere. It is because in 
most parts of the low corona, particularly the strongly-magnetized active regions, the plasma 
/3 (a ratio of gas pressure p to magnetic pressure -B^/(2/io), i.e., /3 = 2fiop/B'^) is extremely 
low (/3 ^ 1) and the plasma velocity v is also low compared to the Alfven speed va {v <^ va), 
which means that the pressure gradient, gravity, and inertial force can be neglected from 
the momentum equation and thus the only-survived Lorentz force must be self-balanced, 
i.e., j X B = (j is the electric current density and B is the magnetic field). This means 
that V X B = aB, where the scalar a is called the force- free parameter. Generally, a varies 
spatially for NLFFF and some popular simplifications include a = for potential field and 
a = constant for linear force-free field. 



The reasons why nonlinear force-free model is superior over other much simpler force- 
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free models for the active regions, i.e., the potential field and the linear force- free field models 



are mainly as follows (Wiegelmann 2008): (1) observation shows that there are significant 
non-potential fields in the active regions, which excludes the potential model; (2) usually 
the force-free parameter a is a very space-dependent function as derived from the measured 
vector magnetogram and also demonstrated by great contradiction of the observed loops 
and linear force-free extrapolations; (3) potential and linear force-free fields are too simple 
to estimate the free magnetic energy and magnetic topology accurately. On the other hand, 
one may wonder why the more realistic model, the MHD model (e.g., (Wu et al. 2006 
2009)), is less commonly used than the NLFFF model. The reason is twofold. Firstly, 



there is a lack of observed information of gas parameters such as the surface plasma density 



and velocity, which are critical boundary conditions for the full MHD simulations (Abbett 



& Fisher 2003 Abbett et al. 2004 Welsch et al. 2004; Wang et al. 2011); secondly the 



numerical realization of full MHD simulation is greatly limited by the present computational 
capability. For instance, the Courant-Friedrichs-Lewy (CFL) condition puts rather severe 
restrictions on the size of the time step for most explicit schemes because densities in the 
corona are very low while magnetic field strengths in active regions can be quite high (~kG) 



which inescapablely results in a extremely high Alfven speed (Abbett & Fisher 2003). The 
computational limitation of full MHD arises especially when applying to very high-resolution 
and large-field-view magnetograms currently available. 

A variety of numerical codes with different methods have been proposed to implement 
the NLFFF extrapolations up to the present. The underlying methods of these codes can 



be classified into six types including (1) the Grad- Rubin method (Grad & Rubin 1958 
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2008 ) for a comprehensive review of most of these methods. In addition to the difference in 
methods, the specific realizations (i.e., the codes) differ significantly in many other aspects 
from software to hardware, e.g., the mesh configuration, the numerical scheme and boundary 
conditions, the language of the code (i.e., IDL, C or Fortran), the hardware architecture and 
the degree of parallelization. As a consequence, these codes perform very differently from 



each other with the computational speed and extrapolation accuracy. Schrijver et al. (2006) 



and Metcalf et al. (2008) have carried out detailed comparisons of some representative codes 
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using respectively the semi-analytic Low & Lou's force- free solutions (Low & Lou 1990) and 



a Sun-like test case constructed by van Ballegooijen et al. (2007). They show that although 



all the tested codes can achieve the reference solutions qualitatively, the differences are 
considerable by quantity. It is pointed out by their analysis that the way of implementation 
of the method plays the same important role as its underlying approach for causing such 



differences. In particular, they found that the optimization method coded by Wiegelmann 



(2004) is the fastest converging and best-performing algorithm. 



In our previous work (Jiang et al. 2011), we have used a fuU-MHD relaxation method 



for reconstructing the corona field basing on our CESE-MHD code (Feng et al. 2007; Jiang 



et al. 2010). We included both the gravity and gas pressure of plasma in the model for a 



more realistic emulation of the low corona. The relaxation is solely dependent on a relatively 
small viscosity term V ■ (z/pVv) (p, v are respectively the plasma density and velocity and 
V is the viscosity) as done by ( [Mikic fc McClymont 1994). It is demonstrated that the 
MHD relaxation method combined with the established CESE-MHD code can gain many 
advantages over other approaches, such as the simplicity of the implementation, the high 
accuracy of the computation, and the efficiency of the highly parallelized code. By the Low & 
Lou's force-free benchmark, it is also found that our implementation reconstructed a result 



comparable with the best one by Wiegelmann (2004) reported in (Schrijver et al. 2006) 



which encourages us for the further work of our method for more realistic applications. As 
another fact, it also proved that the force- free model is a good assumption since we directly 
start from the full MHD and finally reach a state very force-free. Recently in the experiment 
with realistic magnetogram, however, we found that the system is prone to produce very 
large velocity (of several va) when performed on magnetograms with rather large gradient, 
which is ordinary in realistic photospheric data. It is because the discretization errors of 
the large-gradient regions can cause large Lorentz forces whereas the used viscosity v is too 
small to effectively control the motion driven by these forces. This velocity thus severely 
restricts the timestep and further decreases the relaxation speed of the whole system. On 
the other hand, increase of the viscosity may be useful for restricting the plasma velocity, 
but it also significantly restricts the timestep because the CFL condition says 



At < 0.5 



where Ax and At are the mesh spacing in space and time. Too small timestep is in particular 
unfavorable for the CESE scheme, which can produce excessive numerical diffusion and 
lose the accuracy (Chang 2002 Feng et al. 2010). Although an implicit dealing of this 



viscosity term can remarkablely remedy this problem, the price is a big complication of 
the numerical scheme and parallelization. Moreover, for the case of force-free extrapolation 
in which only the magnetic field is solved, it obviously gives no payoff if considering the 
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computational restrictions of the full MHD model, e.g., the slow evolution of the weak 
field and the additional computational resources consumed by solving plasma density and 
pressure. 

In this paper, we propose to present a new implementation of the MHD relaxation 
method for NLFFF extrapolation to avoid the above shortcomings of our previous method. 



We now adopt the magneto-frictional approach as used by (Roumeliotis 1996 Valori et al. 



2007), which explicitly introduce a frictional-like force F = — z/v to the momentum equation. 



By adjusting the frictional parameter z/, one can control the relaxation of the system more 
efficiently than using the viscosity. Different from the convectional form of the magneto- 



frictional equation as used in (Roumeliotis 1996 van Ballegooijen et al. 2000 Valori et al. 



2007) which can not be solved by many modern CFD (computational fluid dynamics) or 



MHD solvers designed for the standard PDE (partial differential equation) system like 

dV dF dG (9H _ g 

dt dx dy dz 



(2) 



we use a form with the time-dependent term of momentum reserved. It will be shown that 
by such modification, the equation system can be still written in standard conservation form 
with source terms, for which the CESE-MHD method is just designed. This paper also 
focuses on a comprehensive examination of the implementation, by applying the code to 



extrapolations of the semi-analytic force-free solutions adopted by Schrijver et al. (2006) 



(hereafter referred to as Paper I) and the more stringent solar-like test used by Metcalf 



et al. (2008) (hereafter referred to as Paper II). All these tests will be carried out with 
the same conditions as much as possible (i.e., the same mesh resolution, the same initial 
potential field, the same artificial boundary conditions) as in the above two papers for a 
rigid assessment and comparison with reported results there. We will show that by this new 



implementation, we successfully improved our method over the previous work of (Jiang et al. 



2011). Quantitative comparisons of the results will demonstrate that our performance of the 



extrapolation accuracy basically reaches the same level of the present best-performing code 



by Wiegelmann (2004) even for the rather stringent test cases. 



The remainder of the paper is organized as follows. In Section [2| we describe the model 
equations and the numerical implementation. In Section |3] we give a briefiy review of the 
benchmark models used for testing the code. The metrics that is used to evaluate the results 
of extrapolations are given in Section |4] and the extrapolation results and comparisons are 
reported and discussed in Section [5j Finally, we offer concluding remarks and some outlooks 
in Section |6l 
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2. The Method 

2.1. The Magneto-Frictional Equations 

In the magneto-frictional method, an artificial frictional force is introduced to the MHD 
momentum balance equation 

Z^v d(pv) 

^ dt + ^ ■ ^'^^^^ = -Vp + pg + J X B - z/v (3) 

where current J = V x B and u is the frictional coefficient. In situations for seeking a force- 
free field, the plasma pressure and gravity can be neglected, which leads to the following 
zero-beta equation 

p = J X B - z/v. (4) 

^ Dt ^ ^ 

By further discarding the inertial term, (i.e., D\ / Dt = 0), it finally gives the usually adopted 



form of the magneto-frictional method (van Ballegooijen et al. 2000 Valori et al. 2007) 



z/v = J X B. (5) 

This is simply a balance between the Lorentz force and the friction term and thus the 
the velocity can be explicitly obtained in terms of the magnetic field. This velocity from 
Equation ^ can then be input to the magnetic induction equation 

— = V X (v X B) (6) 

which drives the evolution of the magnetic field. Note that in such simplification, the only 
equation that needs to be solved is the induction equation and many simple finite-difference 
method can be used to solve it, as long as the frictional coefficient is large enough to suppress 
the potential numerical instability. 

In this paper, we do not use the conventional form of the magneto-frictional equation. 
For convenience of utilizing the existing CESE solver, we partially reserve the inertial term in 
Equation (|4]). Specifically, the time-dependent form of the momentum equation are retained 
as follows: 

^ = (VxB)xB-z.pv, p=|Bp + po, (7) 

where only term V ■ (pvv) are omitted from Equation (|4]). Here the density p is set as for an 
nearly uniform Alfven speed to roughly equalize the speed of evolution of the whole field, and 
a small value po, e.g., po = 0.01 is necessary to deal with very weak field associated with the 
magnetic null. This form is also different from the zero-beta model since the time-variation 
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of momentum p = pv is only induced by the Lorentz force and the frictional term locally 
without affected by the neighboring plasma. This benefits in the context of using a rather 
non-uniform density as p oc -B^. 



For the induction equation, we use 

dB 

~dt 



V X (v X B) - vV ■ B + V(/iV ■ B). 



The terms — vV • B and V(/iV ■ B) added on to the induction equation are both aimed for 
control the numerical error of V ■ B. The first term — vV ■ B is derived from the Powell's 



eight-wave MHD model (Powell et al. 1999) and the second term is a diffusive control of 



V ■ B (Marder 1987; Dedner et al. 2002) with diffusive coefficient /i. The effect of these 



control terms can be explicitly seen by take divergence of the induction equation: 

dpn 



dt 



-V ■ (vp„) + V (yUp„ 



Pr. 



VB. 



(9) 



Equation ^ shows that the numerical magnetic monopoles pm, once arise (either because of 
the numerical error or from the boundary conditions), can not accumulate locally. Instead, 
they are convected effectively with the velocity of the plasma v, and meanwhile is diffused 
among the computational volume with speed of p. 

Another modification is made by utilizing the so-called magnetic field splitting form of 



the MHD equation originated by Tanaka (1994). By dividing the full magnetic field B into 
two parts (B = Bo -|- Bi), a embedded constant field Bq and a deviation Bi, accuracy can 
be gained by solving only the deviation. The magnetic splitting form is usually used for the 
global simulation of the solar wind or its interaction with a magnetized planet such as earth 
( |Tanaka|[l994{ [Nakamizo et~aL][2009j [Feng et al.|[2010| , since in these cases a strong 'intrinsic' 
potential magnetic field is present. In the case of solving a force-free field, a potential field 
that matches the normal component of the magnetogram can be regarded as Bq, which is 
only induced by the current system below the bottom (i.e., the photosphere). While the 
deviation Bi can be seen as the field only induced by the currents in the extrapolation 
volume (above the photosphere). Then the magnetic splitting form of magneto-frictional 
method for solving NLFFF reads as in a complete system 



dpv 

~dt 



(V X Bi) X B - z/pv. 



(9Bi 



V X (v X B) + V(/iV ■ Bi 
9Bo 



vV Bi, 



dt 



= 0, V X Bo = 
p= |B|2 + po, 



O,VBo = 0, 
B = Bo + Bi. 



(10) 
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A notable advantage of using the above equations is that we can totally avoid the random 
numerical currents and divergences remaining in the initial potential field that is computed 
by the Green's function method or other numerical realization. It is commonly made that 
in the extrapolation box the currents are concentrated in the interior of the volume while 
the upper and the surrounding region are dominated by the relatively weak potential field 
(Schrijver et al. 2008 DeRosa et al. 2009). Thus the splitting form can retain the accuracy 
of this field. Other merits of using the splitting equations will be seen in the implementation 
of a multigrid-type optimization (Section 2.2). 



2.2. Numerical Implementation 



The above equation system (10) can be written in a general conservation form with 
source terms as follows 



dV OF dG dU 

dt dx dy dz 



OF, dG, dU, 



dx 



dy 



dz 



with U = (pv,Bi,Bo) and other terms are given in Appendix. We then input this model 
equations to the CESE code, which is designed for any equations that can be written in 
the above standard form. The CESE method deals with the 3D governing equations in a 
substantially different way unlike the traditional numerical methods (e.g., the finite-difference 
or finite- volume schemes). The key principle also a conceptual leap of the CESE method is to 
treat space and time unitedly as one entity. By introducing the conservation elements (CEs) 
and solution elements (SEs) as the vehicles for calculating space-time fiux, the CESE method 
can enforce the conservation laws both locally and globally in their natural space-time unity 
form. Compared to many other numerical schemes, the CESE method can achieve higher 
accuracy with the same mesh resolution meanwhile provides simple mathematics and coding 
such as free of any type of Riemann solver or eigendecomposition. Thus it can benefit for 



the non- hyperbolic system like the present form of magneto- frictional model (10). For more 



detailed descriptions of the CESE method for MHD simulations, please see (Feng et al. 2006 



Zhang et al. 


2006; 


Feng et al. 


2007 


Jiang et al. 


2010, 


2012 



The initial and boundary conditions are given as usually as done in the MHD relaxation 



method for NLFFF extrapolation (Roumeliotis 1996 Valori et al. 2007; Jiang et al. 2011). 



The initial magnetic field is supplied with the potential field Bpot computed using the LoS 
magnetogram. In the magnetic splitting form, it is simply set by Bq = Bpot and Bi = 0. 
The system is started from a static state (v = 0) and driven by inputting vector-magnetic 
information on the bottom boundary. Specifically, at the bottom boundary, the magnetic 



field Bi is changed linearly from the initial value to the final value B^ 



Bpot (B^ 



IS 
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the vector magnetogram) in several Alfven time t^. In such process, the Lorentz forces are 
continuously injected from the bottom to drive the system away from the initial potential 
field. After then the bottom boundary is fixed for the system to relax to a new equilibrium. 
During the whole evolution, the lateral and top faces are fixed as Bi = and the velocity 
of all boundaries is unchanged as v = 0. 

The time step At is restricted by the CFL condition as 

Ax , , 

At = 0.5 (12) 

'VA + '^max 

where the Alfven speed va = and f max is the maximum velocity of the whole computational 
domain. In this context, an arbitrary choice of the frictional coefficient {v > 0) can be 
workable because the numerical instability is prohibited by the CFL condition. But choosing 
a proper u is particularly important since it controls the relaxation speed of the system. A 
simple half-discretizing of the momentum equation ([T]) gives 



n+l _ 

^ ^ JxB-z/p"+i (13) 



At 

where n denotes the time level (note that the source terms, e.g., the friction, are treated 
imphcitly in the CESE method), and thus 



Equation ( 14 ) shows that the effect of the friction is simply to reduce the momentum every 
time-step by a factor of 1 + z/At. A too small z/ is prone to lead to a too large velocity (> va), 
which may distort the field line excessively. On the other hand, a too strong friction will 
suppress the velocity to a very low value that makes the system too hard to to be driven. 
For a compromise we set u = AcAt/Ax"^ which gives the factor 

where c ~ 1 is variable for optimizing the relaxation. In this form, the friction is adaptively 
optimized for both the driving and relaxing processes according to the maximum velocity: 
in the driving process, the velocity is relatively large which thus reduces the friction for fast 
evolution away from the initial field; in the relaxing process, the velocity becomes smaller 
which will increase the friction for fast relaxation to equilibrium. Similar setting of u is also 
done by ( van Ballegooijen et al.|[2000 Valori et al. 2007). Finally for the diffusive coefficient 



of V-B, we set n = 0.4Aa:^/At to maximize the diffusive effect without introducing numerical 
instability. 
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One great challenge of the NLFFF reconstructions is the limitation of computational 
resources, especially for the extrapolation of currently available high-resolution and large- 
field-view magnetograms, thus a parallel computation is generally necessary to be used. Our 
method is parallelized by the AMR-CESE code (Jiang et al. 2010[ Jiang et al. 2011) which 
is a combination of the CESE code within the PARAMESH toolkit (a open-source Fortran 



package for implementation of the parallel-AMR technique on existing code (MacNeice et al. 



2000)), and is performed on a share-memory parallel cluster. Also for a large magnetogram. 



a multigrid-like strategy is recommended to be used to both accelerate the computation 
and improve the quality of extrapolation (Metcalf et al. 2008). We compute the solution 
serially on a number of grids with the resolution ratio of two, and input the results of 
the coarser resolution to initialize the next finer resolution. It should be noted that such 
method is not standard multigrid, since it does not incorporate different grids simultaneously 
and iterate back and forth between coarser and finer grids, but computes the solution of 
different grids only once with the coarser solution used to initialize the finer grid. The 
main advantage by doing this is to give a better (than potential field) starting equilibrium 
on the full resolution grid. Particularly standard node-centered full-weighting restriction 
and prolongation operators of multigrid method are used to transfer data between different 
resolutions. By these operators, any boundary values are interpolated using data on the 
same boundary face and the total flux of the magnetogram is conserved between different 
grids. It is worthwhile noting that when using the coarser results to initialize the finer grid, 
the magnetic splitting form can gain accuracy since only the deviation field Bi is needed 
to interpolate, while the intrinsic potential field Bq is always reset by the value on the 
final resolution or the Green's function method. However, the multigrid algorithm presents 
demerit because every prolongation (by interpolation) will introduce new errors of divergence 
and current to Bi which can be 'felt' by the CESE method. Such problem is also faced in 
many AMR simulations due to the mesh refinement (Toth & Roe 2002). We will discuss its 
effects in Section |5] by comparing the results with and without the multigrid algorithm. 



3. Benchmark Models 



3.1. Low and Lou Force Free Field 



The NLFFF model derived by Low & Lou ( 1990 ) has served as standard benchmark for 



many extrapolation codes (Wheatland et al. 2000 Amari et al. 2006; Schrijver et al. 2006 



Valori et al.|[2007l |He fc Wang|[2008| [Jiang et al.|[20ll| ). The fields of this model are basically 
axially symmetric and can be represented by a second-order ordinary differential equation 



- 11 - 



derived in spherical coordinates 

(1 - cos^ e)-^^^ + n{n + 1)P + a2i±^pi+2/« = Q (16) 
where n and a are constants. Then the magnetic field is given by 

r"' sm ay r sm u or r sm 



where A = P(cos0)/r" and Q = aA^^^I"^. The solution P of Equation (16) is uniquely 



determined by two eigenvalues, n and its number of nodes m ( Low fc Lou||199d Amari et al. 



2006). By arbitrarily positioning a plane in the space of the analytical fields, one obtains 
a different kind of test case, in which the plane represents the bottom-boundary condition 
for extrapolation of the overlaying fields. In this way the fields sliced by the plane show no 
more symmetry and thus benefit for a general testing of extrapolation. The position of the 
plane is characterized by two additional parameters, / and $. Here we choose two particular 
solutions characterized by the parameters (n, m, /, $), which are respectively given by n = 1, 
m = 1, I = 0.3 and $ = 7r/4 (referred to as CASE LLl), and n = 3, m = 1, I = 0.3 and 
$ = 47r/5 (CASE LL2). For both cases, the computational domain is x,?/ G [— 1,+1] and 
z G [0, 2] and discretized by uniform grid of 64 x 64 x 64 (same as Paper I). The same test 
solutions are also used by works in the above references where more analyses of these fields 
can be found. The vector magnetograms for both cases at z = are shown in Figure [T] 
and their 3D field lines are shown in panels (a) of Figures |4] and [7j Basically, non-potential 
fields occupy more volume in CASE LLl than CASE LL2 and CASE LL2 is 'more nonlinear' 
with a larger a and stronger fields more concentrated near the center of the model volume. 
Since the solutions show rather smooth (small gradients) and relatively simple magnetic 
structures with their topologies roughly consistent with those of the potential fields based 
on the same LoS magnetograms, these tests can be regarded as preliminary tests for any new- 
developed NLFFF extrapolation methods before facing to more stringent cases or realistic 
magnetograms . 



3.2. The van Ballegooijen Reference Model 

The van Ballegooijen reference model is adopted from Paper II for a more stringent and 
realistic testing of our code. By this reference model, it is possible to mimic the analysis of 
real observational data while still knowing the properties of the field to be modeled. This 
model field is constructed by initially inserting an S-shaped fiux bundle into a potential field 
associated with active region AR 10814 (see panel (a) of Figure [3]), and then relaxing the 
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unbalanced system to a near force-free state using van Ballegooijen's magneto-frictional code 
in spherical geometry ( van Ballegooijen et al.||2000 van Ballegooijen||2004 ). Furthermore, an 



upward force was apphed to the field at the lower boundary during the relaxation process to 
mimic the effect of magnetic buoyancy in the photoshpere, and thus achieve more realistic 
magnetic fields between photospheric and chromospheric heights in the model. Finally, by 
coordinate transformation and interpolation from the original spherical geometry, magnetic 
fields in a Cartesian box of 320 x 320 x 258 pixels centered on the active region is extracted 
as the final reference model. The final near-force-free field are drawn in panel (b) of Figure [3| 
which contains several interesting topological features including a coronal null and its associ- 
ated separatrix surface, and the S-shaped flux bundle surrounded by a quasi-separatrix layer 
(please see more details in Paper II). Compared with the initial potential field, the magnetic 
topology near the bottom is significantly modified by the low-lying flux rope, which chal- 
lenges the extrapolation much more than the Low and Lou cases. Because of extra force 
presented at the bottom, this reference model can be used for tests of extrapolations from ei- 
ther the 'chromospheric' or the 'photospheric' magnetograms, by providing the NLFFF code 
with data ok, z = z^oi z = Zq [z is the height in the model, e.g., z^ is the base of the reference 
model). For the 'chromospheric' case, the boundary data used is largely force-free which is 
consistent with the extrapolation method, while the 'photospheric' case is more forced and 
thus represents a more realistic magnetogram of observation. It should be stressed that 
the van Ballegooijen reference model is not strictly force-free in the whole model box, even 
above the chromosphere, due to the implementation of the magneto-frictional method and 
some other numerical errors. It is demonstrated that in the model the residual forces of at 
least 5% of magnetic-pressure force are present up to height of Z30, which are consistent with 
what is known of forces on the Sun, making the model a realistic, solar-like test case for the 
extrapolation codes. 

As done in Paper II, we will test our code by both the chromospheric and photospheric 
cases. For the photospheric case which is inconsistent with force-free assumption, we only 



examine the code with the photospheric magnetogram preprocessed by method of Wiegel- 



mann & Neukirch ( 2006 ) . The magnetograms of both test cases are plotted in Figure [2| 



which show a significant shearing along the polarity inversion line (PIL). 



4. Metrics 



For a detailed analysis of the extrapolation fields, a suite of metrics introduced in (Schri- 
jver et al. 2006) are computed. These metrics compare either local characteristics including 
vector magnitudes and directions at each point or the global energy content. They are 
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respectively the vector correlation C^^c 

i i 

the metric Ccs based on the Cauchy-Schwarz inequality 



BjIbJ 



(19) 



the normalized and mean vector error E'^, E'^ 

E^ = Y, |bi - B,|/ ^ |B,|; K = 1 - ^n, 

i i 

1 \ ^ |Bi — hi 



Er, 



IB. 



-■ E' = 1 - E 



(20) 
(21) 



where Bj and bj denote the input field (the Low & Lou's solution or the van Ballegooijen 
reference model in this paper) and the extrapolated field, respectively, i denotes the indices 
of the grid points and M is the total number of grid points involved. As can be seen, an exact 
extrapolation will have all the metrics equal to unity in such definitions, and the closer to 
unity means the better extrapolation and vice versa. Detailed descriptions for these metrics 



can be found in (Amari et al. 2006; Schrijver et al. 2006; Valori et al. 2007). Another very 



important parameter for comparing the extrapolation is the free energy of magnetic field. 
It is measured by the ratio of extrapolated energy to the potential energy using the same 
magnetogram, 

EJB.p 



E/E, 



pot 



(Bpot)ip 



(22) 



It is common to measure the force-freeness and divergence-freeness of the extrapolation 



using the current- weighted sine metric CWsin and divergence metric (Metcalf et al. 



2008 


Schrijver et al. 


2008 


DeRosa et al. 


2009 


Wheatland et al. 


(2000 


) as 



CWsin 



Ei l^iWi 



\3i X Bi 



and 



m) 



EJJ, 

£^ (V-B 



Jill Bi 



M ^ 6|B,|/Aa; 



(23) 



(24) 



where Ax is the grid spacing. Both of the metrics are normalized with the former focused on 
the directional deviation between the currents and the field lines and the latter on the relative 
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value of residual divergence. These two metrics are equal to zero for an exact force-free field, 
and hence the smaller these metrics are, the better the extrapolation is. 

In addition to the above metrics, we also introduce another pair of metrics to evaluate 
the degree of convergence towards the divergence-free and force-free state. For the first one, 
we note that a nonzero V • B (i.e., the magnetic monopole) introduces to the system an 



unphysical force F = BV • B parallel to the field line (Dellar 2001). To evaluate the effect 
of this unphysical force to the numerical computation, the metric -Ev b is defined as the 
average ratio of this force to the magnetic-pressure force 

_ 1 |B.(V-B).| 
^v-B-^l^|^(|3|2/2)^|. (25) 

Similarly, the second metric -Evxb measures the effect of the residual Lorentz force in the 
same way 



1_ sr-^ \3i X Bi\ 

U 2^ i^nRis/ov.r ' 



£^VxB M ^ |V(|B|V2) 

Unlike the metrics of CWsin and which mainly characterize the geometric properties 
of the field, these two metrics directly measure the physical action of the residual divergence 
and Lorentz force on the system in the actual numerical computation. This is important for 
checking of the NLFFF solution if it is used to initialize any MHD simulations. 

For the all metrics above, the first four are more rigid since they are involved without 
any type of derivatives, while the other metrics may be unreliable for comparison with results 
from different papers due to the specific numerical realization of the derivatives (for example, 
different orders of numerical differentiation or different configurations of computational grid, 
e.g., cell-centered or staggered). In the present work, the second-order central difference 
is used for evaluating all the derivatives associated with the divergence, curl, and gradient 
operators, although the spatial derivatives can be directly obtained from the CESE method. 



5. Results 



In this section, we present the results of extrapolation for the benchmark models. The 
results are also compared with some results reported in Paper I, II and (Valori et al. 2007). 
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5.1. Low and Lou's Force Free Field 

5.1.1. CASE LLl 

Results for CASE LLl are given in Figure |4], Figure [5} Table [T] and Table [2] In Figure |4| 
we show the same selected field lines in 3D view for the extrapolation results and the reference 
solution, which are traced from foot-points evenly rooted at the lower boundary. In the 
central region of core fields (i.e., x,y E [—0.5, +0.5] and z G [0,1], enlarged in the bottom 
row of the figures), the MHD result is highly in agreement with Low & Lou's solution, as can 
be seen from the high similarity of most of the field lines. Such agreement is demonstrated 
quantitatively by the metrics in Table [Tj The first three metrics are very close to 1 with 
error of < 5% and even the most sensitive metric E'^ has an error below 10%. In Table |l| 



we also compare with the best result by Wiegelmann's code (Wiegelmann 2004) reported in 



Paper I and extrapolation by Valori et al. (2007). Our result for the central region, although 



only specified the lower boundary, still reaches the level of the best extrapolation that using 
information of Low & Lou's solution on all six boundaries. This may be due to the fact that 
this test case is close to potential field hence a fixed side and top boundary conditions can 
rarely impact the central extrapolation. The infiuence of the boundary conditions is more 
explicitly shown by comparison of the metrics for the entire domain. Now the Wiegelmann's 
extrapolation still performs perfectly with all four metrics extremely close to unity, while our 
result and the Valori's perform less exactly, but still satisfactory. Like other results in the 
Table, our extrapolation also recovered the energy content very precisely, especially for the 
central region; furthermore the metrics evaluating the force-freeness and divergence-freeness 
are rather small and close to the reference values which is caused by discretization error. All 
these shows that extrapolation of very high accuracy can be achieved by our implementation, 
at least for the present, relatively easy test case. 



In the Valori's implementation of the magneto-frictional method (Valori et al. 2007), 
a fourth-order numerical scheme is used with many layers of ghost mesh and a high-order 
polynomial extrapolation of the fields is adopted on the side and top boundaries. However by 
comparing our result with Valori's, it is interesting to note that our implementation performs 
better although our numerical scheme is second-order method without any ghost layers and 
the boundaries are simply fixed. In this test, the boundary effect can be neglected as said 
(while in the following test of CASE LL2, we will see the effect of different treatments of 
these boundaries). Then we concluded that the better performance is due to the merit of 
using a better solver, i.e., CESE method combined with the magnetic splitting algorithm 
which can gain additional accuracy. 

We finally give a study of the convergence of the extrapolation. In Figure |6] we shows 
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the history of the system relaxing to the final force-free equilibrium, including the residual 
of temporal evolution of the magnetic field 



res 



(B) 



n— 1^ 
5 , 



=x,y,z 



Em. 



(27) 



(where n denotes the iteration steps), the evolution of the velocity and the metrics. The 
system converged very fast from a initial residual of 10^^ to value under 10^^ with time 
of lOOr^ (about 5000 iteration steps, see panel (a) of Figure [6]). The evolution of the 
plasma velocity indicates that a static equilibrium is reached as expected with a rather 
small residual velocity ~ 0.01 which is only on the order of the numerical error OlAx"^) of 
the CESE solver. All the metrics plotted in the figure converged after AOta (about 2000 
iteration steps), when the residual is on the order of 10~^. Note that the metric /j of V ■ B, 
like the plasma velocity, first climbs to a relatively high level (~ 0.01, see panel (d)) and 
then drops to the level of discretization errors. In principle the divergence-free constraint 
of B should be fulfilled throughout the evolution, at least close to level of discretization 
error. However, a ideally dissipationless induction equation Equation ^ with divergence- 
free constraint can preserve the magnetic connectivity, which makes the topology of the 



magnetic field unchangeable (Wiegelmann 2008) unless a finite resistivity is included for 



allowing reconnection and changing of the magnetic topology (Roumeliotis 1996). In the 



present implementation in which no resistivity is included in the induction equation, a break 
of V ■ B constraint in the initial evolution process (indicated by the climb of metric /«) can 
thus produce a change in the magnetic topology (also note that a numerical diffusion can 
help to topology adjustment in addition). 



5.1.2. CASE LL2 

Now we present the result of the second test CASE LL2 which is more difficult than 
the first one. Compared with CASE LLl, this case is more non-potential and nonlinear with 
relatively larger gradient of fields. Although considering of this, our extrapolation still gives 
satisfactory result that is as good as the best result in Paper I (see the first four metrics 
shown and compared in Table |3]) . It is also quite encouraging that our result for the central 
region is close to Valori's which is computed with fourth-order numerical scheme. Figure [7] 
and [8] demonstrate qualitatively the consistence with the Low and Lou's reference solution. 
The energy contents are very well reproduced for both the central region and the entire 
domain. 

By comparing the metrics of the entire domain, we find that our result scores worse than 
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Model acc Ccs E'^ E/E^o, CWsin (|/,|) 

For the central region 
Low 

Our result 

Wiegelmann* 
Valori** 
Potential 



For the entire domain 
















Low 


1.000 


1.000 


1.000 


1.000 


1.294 


0.014 


0.56 X 10-^ 


Our result 


0.998 


0.955 


0.873 


0.662 


1.282 


0.060 


1.31 X 10"^ 


Wiegelmann* 


1.00 


1.00 


0.98 


0.98 


1.31 


0.070 




Valori** 


0.994 


0.86 


0.80 


0.51 


1.28 


0.019 




Potential 


0.852 


0.824 


0.446 


0.353 


1.000 







1.000 


1.000 


1.000 


1.000 


1.242 


1.000 


0.997 


0.964 


0.912 


1.241 


1.00 


1.00 


0.97 


0.96 


1.26 


0.999 


0.99 


0.95 


0.87 


1.23 


0.858 


0.869 


0.498 


0.443 


1.000 



0.014 
0.015 

0.009 



0.94 X 10"^ 
1.67 X 10-^ 



Table 1: CASE LLl: Metrics for the central and entire regions. The superscript * denotes 
the reported results in Table I and Table II of Paper I, and ** denotes the reported results 
bylValori et all (|2007l). 



Model 






For the central region 






Low 


0.007 


0.005 


Our result 


0.010 


0.008 


For the entire region 






Low 


0.003 


0.002 


Our result 


0.024 


0.005 



Table 2: CASE LLl: Metrics of Evxb and Ey.B 
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the Valori's, especially shown by E'^. The reason for this is twofold. Firstly, a high-order 
scheme can characterize the large gradient much more accurately than 2nd-order scheme and 



thus the fourth-order scheme by Valori et al. (2007) shows its advantages for this test case 



with larger gradient than CASE LLl (the high-order accuracy can be also achieved by mesh 
refinement other than improving the order of the numerical scheme. As demonstrated in our 
previous work ( Jiang et al.|20H ), a refined grid of 128 x 128 x 100 with the CESE scheme gave 
much more accurate result with the most sensitive metric even reaching 0.166). Secondly, our 
implementation simply fixed the artificial boundaries (i.e., the lateral and top faces), which 
obviously make the system over-determined. In the present test case the final field is very 
non-potential, hence the boundary values deviate from the initial conditions significantly. 
A fixed boundary condition will tie the field lines that pass through the boundary, which 
would otherwise freely cross the boundary since this boundary is a non-existing interface in 
the realistic corona. This line-tied condition thus hinders the system from relaxing to a true 
force-free state. The metric of CWsin demonstrates clearly that the field is less force-free 
than the Valori's result. The boundary effect can be seen visually from some field lines close 
to the lateral boundaries shown in the figures (e.g., see the distortion of field lines at the 
lower left coroners of panel (b) in Figure [T] and |8]) . Also caused by the boundary effect, the 
field for the entire domain (CWsin = 0.060) is thus less force-free than the central region 
(CWsin = 0.047). In a full MHD simulation, this boundary effect can be minimized by using 
a so-called non-reflecting boundary conditions based on characteristic decomposition of the 



(Wu et al. 


2001 


2006 


Hayashi 


2005 


Feng et al. 


2010 


Jiang et al. 


2011 



However for the present form of magneto-frictional equation, the characteristic method is 
no more valid because of the eigen degeneration of the non-hyperbolic system. In such 
situation, a natural choice of modeling the non-reflecting boundary is to use linear or high- 



order extrapolation, just as done by Valori et al. (2007). Fortunately, field configuration 



like this Low & Lou case with the entire volume very non-potential is not usually found in 
observed magnetograms, and by choosing a model box significant larger than the core non- 
potential region, a simply fixed boundary values of the initial potential field is still sufficient 
for most extrapolations, for example, test cases in the next section. 

Figure |9] shows that the convergence speed is even faster than case LLl. The residual 
of the magnetic field reaches the order of 10"'' with only ~ 3200 iterations at QQta- The 
evolution speed of the magnetic field is also demonstrated by magnitude of the plasma 
velocity, which is larger than that of case LLl. At about SOta, all the metrics converged 
with the corresponding residual of magnetic field on the order of 10~^. This and the preceding 
cases, shows that generally the iteration can be stopped when res"(B) < 10~^ since after 
that temporal change of the magnetic field can be neglected actually. 



- 19 - 



Model 


C 


Ccs 


K 




E/ Epot 


CWsin 


m) 




For the central region 


















Low 


1.000 


1.000 


1.000 


1.000 


1.099 


0.036 


3.14 X 10" 


-4 


Our result 


0.999 


0.933 


0.938 


0.636 


1.114 


0.047 


7.35 X 10- 


-4 


Wiegelmann* 


1.00 


0.91 


0.92 


0.66 


1.14 








Valori" 


0.999 


0.95 


0.96 


0.75 


1.11 


0.015 






Potential 


0.923 


0.661 


0.572 


0.299 


1.000 








For the entire domain 


















Low 


1.000 


1.000 


1.000 


1.000 


1.100 


0.035 


2.17 X 10" 


-4 


Our result 


0.999 


0.574 


0.852 


-1.060 


1.115 


0.060 


6.35 X 10" 


-4 


Wiegelmann* 


1.00 


0.57 


0.86 


-0.25 


1.14 








Valori** 


0.999 


0.64 


0.88 


-0.01 


1.11 


0.027 






Potential 


0.921 


0.346 


0.465 


-0.639 


1.000 









Table 3: Same as Table [T] but for LL CASE2. The superscript * denotes the reported results 



in Table I and Table II of Paper I, and ** denotes the reported results by Valori et al. (2007). 



Model 


-E'vxB 


-E'v-B 


For the central region 






Low 


0.011 


0.008 


Our result 


0.026 


0.024 


For the entire region 






Low 


0.004 


0.003 


Our result 


0.025 


0.013 



Table 4: CASE LL2: Metrics of EyxB and ^v b 
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5.2. The van Ballegooijen Reference Field 

We first perform tlie extrapolation of the chromospheric case, which provides a largely 



force-free magnetogram without any preprocessing. Figure 10 shows the 3D field lines of 
extrapolation and the reference model for a side-by-side comparison. The field lines shown 
are traced from foot-points equally spaced at the bottom surface. We specially adjusted the 
figures to an approximately same orientation as Paper II on purpose of visual comparison 
with those results by other codes. The same field lines are also projected on the x-y plane in 



Figure 11 As shown from a overall view of the figures, the extrapolation reproduced quite 
well the basic magnetic topology including the low-lying S-shaped field-line bundle, and the 
overlying magnetic arcade that straddles the flux rope. For a confirmation of presence of 
the flux rope in the extrapolation, we select a set of field lines of the S-shaped bundle and 



plot them using different colors in Figure 12 with different perspectives. It clearly shows 



that the flux rope are qualitatively recovered, encouraging us that the code can be used to 
handle such relatively complex test cases. The extrapolated flux rope is weakly twisted and 
its core fields basically lies along the bottom PIL, showing a highly shear with respect to the 
overlying arcades. 

In Table [5] we give a quantitative evaluation and comparison of the result with those 
reported in Paper II. To minimize the (side and top) boundary effects, the metrics are applied 
to the central region of {x^y) G [48,271] x [48,271] with two different heights: z G [2,61] 
for focusing on the low- lying flux rope and z G [2, 225] for including also the surrounding 
potential-like field. This is done in the same way as in Paper II and the results along with 
the first three preferable results in Paper II are presented. Note that some metrics (e.g., the 
Ccs and CWsin) for the reference and potential models differs slightly from those shown in 
Paper II, which may be due to different precisions (i.e., single or double precision) used in 
the numerical computation. 

As shown by Table [sj the extrapolation performs very well with errors only of ~ 2% for 
the first two metrics. These metrics, however, are not sensitive indicator of extrapolation 



accuracy (Valori et al. 2007) and thus all the methods in Paper II give these metrics within 
the same level, even including the potential field. By the more sensitive metrics and 
we find that our result for z G [2,225] is identical with the best result performed by 
Wiegelmann. For the lower height, the deviations of these metrics with the best result 
are smaller than < 2% (note that these deviations can be also partially introduced by the 
different precisions in the numerical computation). This comparison encourages us further of 
the success of our implementation. Our result also gives the energy metrics very close to the 
reference values, showing a good recovering of the free energy content, which is particularly 
important for coronal field extrapolation. Finally the CWsin metric behaves a bit worse 
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than the first five metrics if compared with the Wiegelmann's result, but it still scores better 
than the following results by Wheatland and Valori. 

The last two metrics -Evxb and -Ey-B (shown in Table]?]), both of which are non-zero but 
in a very low level, mean that the reference model is very close to force-free and divergence- 



free but never exact as expected (see Section 5.2). Again, our results present values very 



similar to those of the reference model. Note that although both CWsin and -Evxb nieasure 



the degree of force-free, the former depends more strongly on the high-current regions (Valori 
et al.||2007 ) while the latter is not. Hence CWsin for both heights gives nearly the same value 
of 0.1, while xB gives significant larger value 0.06 for the lower domain [z G [2,61]) than 
value 0.018 for the full domain [z G [2, 225]), which shows that the residual forces are mainly 
presented in the lower region. By these two metrics, we demonstrate that our code can 
minimize the residual forces into a very low values and fulfill the divergence-free condition 
well. 

In Table [5| the result with the multigrid-type algorithm is also presented, which is 
performed on a three-grid sequence (80 x 80 x 64, 160 x 160 x 128 and the full resolution 
320 X 320 X 256). We found that this result (with multigrid algorithm) behaves a little bit 
worse compared with the result without multigrid algorithm. This is due to the additional 
numerical errors of V x B and V ■ B introduced by interpolating field from coarser to finer 
grid, which thus results in a little bit larger values of CWsin and than those without 
multigrid algorithm. Without the multigrid algorithm, the magnetic field splitting form can 
avoid such numerical errors by using a zero Bi initially and thus performs better. However 
in spite of such small disadvantage, the multigrid algorithm is still encouraged to be adopted 
since its reward is a significant reduce of the computing time (e.g., for the present test case, 
we find that using of the multigrid algorithm saves approximately two-thirds of the CPU 
time) . 

Now we present the results for the photospheric case with preprocessed magnetogram. 
The preprocessing procedure can remove the net forces of the photospheric magnetogram 
and hence provide the extrapolation code a more consistent lower boundary condition than 
the raw magnetogram. In Paper II, both cases with the raw and preprocessed photospheric 
magnetograms are tested, and it is found that the preprocessed case gives a significantly 
better result than the raw case in which the fiux rope is not reproduced at all. This demon- 
strated that the preprocessing procedure is necessary for those extrapolation codes. Besides, 
the smoothing involved in the preprocessing also benefits the extrapolation code which is 
based on numerical finite-difference. 



Our result are given in Figure 13 and Table [6] The field lines shows that the overall 
structure including the fiux rope is still recovered qualitatively, but partially. By careful 
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comparison with the reference or the chromospheric cases, the difference is also evident, e.g., 
the bundle of S-shaped lines is much thinner in the present case, which seems that many 
field lines have not reconnected fully to form a entire S. The quantitative comparison also 
demonstrates that the extrapolation quality is a bit worse than chromospheric case, especially 
shown by the CWsin and -EvxB metrics (nearly twice of those of the chromospheric case). 
It means the extrapolated field is farther away from the exactly force-free state than the 
chromospheric case. The free energy is also underestimated much, as the results in Paper II. 
This is because the preprocessing does not recover a force-free magnetogram that is entirely 
consistent with the reference model. Still, it is worthwhile noting that our result again 
reaches the same level of the best in Paper II, and some of the metrics, including the most 
sensitive one E^, perform even better than Wiegelmann's. 



Model 


r' 


Cos 


K 


K 


E/ Epot 


CWsin 


m) 




For ze[2, 225] 


















Reference 


1.000 


1.000 


1.000 


1.000 


1.343 


0.107 


1.02 X 10- 


-4 


Our result 


0.995 


0.978 


0.885 


0.734 


1.337 


0.133 


1.11 X 10- 


-4 


Our result^ 


0.993 


0.975 


0.868 


0.712 


1.357 


0.142 


1.32 X 10- 


-4 


Wiegelmann* 


1.00 


0.99 


0.89 


0.73 


1.34 


0.11 






Wheatland* 


0.95 


0.98 


0.79 


0.70 


1.21 


0.15 






Valori* 


0.98 


0.98 


0.84 


0.71 


1.25 


0.15 






Potential 


0.852 


0.952 


0.687 


0.665 


1.000 








For z e [2, 61] 


















Reference 


1.000 


1.000 


1.000 


1.000 


1.361 


0.103 


1.33 X 10" 


-4 


Our result 


0.995 


0.989 


0.923 


0.883 


1.348 


0.125 


1.58 X 10- 


-4 


Our result^ 


0.993 


0.979 


0.902 


0.825 


1.366 


0.129 


1.74 X 10- 


-4 


Wiegelmann* 


1.00 


0.99 


0.94 


0.89 




0.11 






Wheatland* 


0.95 


0.95 


0.79 


0.76 




0.15 






Valori* 


0.98 


0.96 


0.86 


0.81 




0.15 






Potential 


0.847 


0.901 


0.660 


0.678 


1.000 









Table 5: The van Ballegooijen reference model: metrics for extrapolation of the chromo- 
spheric case. The superscript M denotes the multigrid-type algorithm is used for speeding 
up of the extrapolation and the superscript * denotes the reported results in Table 3 and 
Table 4 of Paper II. 
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Model 


^vec 




E' 

n 


E' 

m 


1 -^pot 


CWsin 


\U* 1/ 




For ^ e [2, 225] 


















Reference 


1.000 


1.000 


1.000 


1.000 


1.531 


0.107 


1.02 X 10^ 


-4 


Our result 


0.970 


0.968 


0.783 


0.679 


1.146 


0.257 


2.29 X 10" 


-4 


Wiegelmann* 


0.98 


0.97 


0.77 


0.65 


1.18 


0.26 






Wheatland* 


0.88 


0.96 


0.69 


0.65 


1.03 


0.11 






Potential 


0.850 


0.945 


0.659 


0.636 


1.000 








For ^ e [2, 61] 


















Reference 


1.000 


1.000 


1.000 


1.000 


1.559 


0.103 


1.33 X 10" 


-4 


Our result 


0.970 


0.980 


0.800 


0.791 


1.149 


0.230 


3.45 X 10- 


-4 


Potential 


0.845 


0.891 


0.629 


0.646 


1.000 









Table 6: The van Ballegooijen reference model: metrics for extrapolation of the photospheric 
case. The superscript * denotes the reported results in Table 3 of Paper II. 



Model EvxB ^VB 

For z e [2, 225] 
Reference 

Our result of the chromospheric case 
Our result of the photospheric case 
For z e [2, 61] 

Reference 0.060 0.010 

Our result of the chromospheric case 0.072 0.014 

Our result of the photospheric case 0.132 0.040 

Table 7: The van Ballegooijen reference model: metrics of -Bvxb and -Bv b- 



0.018 0.004 
0.023 0.005 
0.048 0.015 
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Conclusions 



As a viable way of study magnetic field in the corona, the NLFFF extrapolation also 
needs a considerable effort to be devoted for its numerical realization. In this paper, a new 
numerical implementation of NLFFF extrapolation is presented based on the MHD relaxation 
method and the CESE-MHD code. Our implementation outstands for the following aspects. 

1. The magneto-frictional approach that is designed for speeding the relaxation of 



^Roumeliotis 


1996 


Valori et al. 


2007 



performance CESE scheme on a grid without any type of ghost zone or buffer layer. 

2. The accuracy is further improved by firstly utilizing the magnetic splitting form 



(Tanaka 1994) for NLFFF methods to totally avoid the numerically random errors involved 



with the initial input. 

3. Multi-method control, i.e., the diffusive and convection term, of numerical magnetic 
monopoles is employed for effectively reducing the divergence error. 

4. The vector magnetogram is inputted at the bottom boundary in the way of time- 
linearly modifying the potential field to matching the magnetogram and other artificial 
boundaries are just fixed as the initial potential values, which makes our implementation 
much easier than other MHD relaxation methods (e.g., those by Roumeliotis (1996) and 



Valori et al. (2007)). 



5. The code is highly parallelized with the help of PARAMESH toolkit and performed 
on the share-memory parallel computer. It can be readily realized with the AMR technique 
and applied to very high-resolution magnetogram in the near future. The multigrid-type 
algorithm is also incorporated into the code to speedup the computation as recommended. 

We have examined the capability of the method by several reference solutions of NLFFF 
that can be served as a suite of benchmark tests for any NLFFF extrapolation code. These 
test cases consist of the classic half-analytic force-free fields by Low & Lou ( Low fc Lou|1990 ) 



and the much more stringent and solar-like reference solution by van Ballegooijen (2004). 
The results show that our method are successful and versatile for extrapolations of either the 
relatively simple cases or the rather complex cases which need significant rebuilding of the 
magnetic topology, e.g., the flux rope. We also compute a suite of metrics to quantitatively 
analyze the results, and shows that the solutions are extrapolated with high accuracies which 



are very close to and even surpass the best results by Wiegelmann (2004) from comparison of 
the metrics. This demonstrated that, at least in computation accuracy, our code performs as 
good as the best one of state-of-the-art (the computing time of the code, however, is difficult 
to be compared because the hardwares are different). In addition, we introduced a pair of 
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metrics for assessment of the divergence-freeness and force- freeness of the extrapolation, -Ev b 
and i^vxB; which further demonstrated that our code can fulfill the solenoidal constraint 
well and minimize the Lorentz force to the same level of the reference values. 

The success of our implementation encouraged us that with a good solver, the MHD re- 
laxation approach can also extrapolate the NLFFF as accurately as other good-performance 
algorithms, like the weighting optimization method. This confirms again that the way of 
implementation of the methods plays the same important role as their underlying approach. 
It is especially worthwhile pointing out that, as also noted by ( Wiegelmann||2008 ) , the MHD 
relaxation approach has a great advantage over other methods: any available time-dependent 
MHD code can be adjusted for NLFFF extrapolations, which thus saves the major effort 
that should be made to develop a new code from scratch for a special method. We are 



also inspired by Valori et al. (2007), who show that a higher order scheme can significantly 



advance the extrapolation. In our project, an arbitrary high order CESE scheme is under 
development and is expected to be used for future improvement of our implementation. 



Recently, more critical tests of extrapolation codes have been performed by Schrijver 



et al. (2008) and DeRosa et al. (2009) based on vector magnetograms of AR 10930 and 10953 



from Hinode/SOT and observed coronal loops. It is found that the Grad- Rubin-style current- 



field iteration implemented by Wheatland (2006) surpassed the Wiegelmann (2004)'s code 



that performs best in the benchmark tests, and basically the results by different methods are 
very inconsistent with each other. This shows that the idealized tests are unable to assess 
the code's ability to deal with various uncertainties or errors in the real magnetograms and 
more critical assessment of the code using realistic vector magnetograms is also planned in 
our future work. 

The present extrapolation in Cartesian geometry is often limited to relatively local areas, 
e.g., a single active region without any relationship with others. However, the active regions 
cannot be isolated since they generally interact with neighboring ARs or overlaying large 
scale fields. It is also pointed out that the fields of view in Cartesian box are often too small 



to properly characterize the entire relevant current system (DeRosa et al. 2009). To study 



the connectivity between multi-active regions and extrapolate in a larger field of view, it is 
necessary to take into account the curvature of the Sun's surface by extrapolation in spherical 
geometry partly or even entirely, i.e., including the global corona ( Wiegelmann|2007 Tadesse 
et al.|[20"lla||b ). Moreover, a global NLFFF extrapolation can also avoid any lateral artificial 
boundaries which are inescapable and cause issues in Cartesian codes. We are now on the 
way of developing a global NLFFF code for the new era of routinely observation of the 
global vector magnetogram (which will be opened by SDO/HMI). Recently in a project of 
constructing a data-driven MHD model for the global coronal evolution, we have established 
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the CESE method on a so-called Yin- Yang overlapping grid in spherical geometry (Kageyama 



&: Sato|[2004 ). This implementation, combined with our present NLFFF code, will make the 



realization of a global NLFFF extrapolation very viable if provided with the global vector 
magnetogram. The Yin- Yang grid is composed of two identical component grids that are 
combined in a complemental way to cover an entire spherical surface with partial overlap on 
their boundaries. Each component grid is a low latitude part of the latitude-longitude grid 
without the pole and hence the grid spacing on the sphere surface is quasi-uniform. In this 
way, we can avoid the problem of grid convergence or grid singularity at both poles which 
will otherwise arise if an entire spherical-coordinate grid is used, as Wiegelmann (2007) has 
pointed out. However, up to now, there is no suitable test case used for the global NLFFF 
extrapolation other than the simple axially-symmetric Low & Lou cases. Most recently. 



Contopoulos et al. (2011) give a variety of global near force-free solutions by a force-free 



electrodynamics code, using solely the radial magnetogram. Their solutions, however, are 
not unique to the same radial magnetogram, but depends on the initial conditions and on 
the particular approach to steady-state Anyway, we believe their solutions can be used 
as much more realistic and solar-like tests for global NLFFF codes than the semi-analytic 
solutions. In our future work, we will develop and test a new global NLFFF code by the 
global force-free solutions from (Contopoulos et al. 2011). 
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^They are improving their code to also make use of the vector magnetogram to define the solutions 
uniquely (by private communications) 



A. The specific form of Equation (11) 



Bi ■ Bi/2 — Bf^ — {BqxBix + Bi^^Bqx) + Bq • Bi 
—BixBiy — (BQ^Biy + Bi^Boy) 
~BixBiz — {BqxBiz + Bi^Bqz) 




v^B^ - v^B^ 






—BiyBix — {B^yBix + BiyBox) 
Bi ■ Bi/2 — Bly — {BoyBiy + BiyBoy) + Bo ■ Bi 
— BiyBiz — {BoyBiz + BiyBoz) 



^y-Bx '^xBy 



VyBz - VzBy 
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— BizBix — {BqzBix + BizBqx) 

—BizBiy — {BozBiy + BizBoy) 
Bi • Bi/2 — Bf^ — {BqzBiz + BizBqz) + Bq • Bi 



H = 



VzBx 



VxBz 



VzBy — VyBz 








(Al) 



= (0,0,0, /xV-Bi, 0,0, 0,0, Of ; 
= (0,0,0,0,//V-Bi,0,0,0,0)^; 
H, = (0, 0, 0, 0, 0, AtV • Bi, 0, 0, 0)^; (A2) 



S = {-vpvx, -lypvy, -upv^, ■ Bi, VyV ■ Bi, VzV • Bi, 0, 0, 0) (A3) 
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Fig. 1. — Vector magnetograms of the central region x, ?/ G [—0.5, +0.5] for CASE LLl (left) 
and CASE LL2 (right). The contours represent Bz- The tangential field is shown by the 
vectors with blue color in positive region and red in negative B^ region. 
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Fig. 2. — Vector magnetograms of the central 224^ pixels for chromospheric case (left) and 
preprocessed photosplieric case (right). The contours represent with a saturation level 
of ±200 G. The tangential field is shown by the vectors (plotted at every second grid point) 
with blue color in positive region and red in negative B^ region. 



- 31 - 




Fig. 3. — The van Ballegooijen reference model: magnetic field lines with contour of on 
the bottom surface. The field lines shown are traced from footpoints equally spaced at the 
bottom surface, (a) the initial potential field, (b) the final near-force-free reference model. 
The bottom row enlarges the central region outlined by the small cube in the top row. 
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Fig. 4. — CASE LLl: magnetic field lines with contour of Bz on the bottom surface, (a) 
the Low & Lou's solution, (b) the extrapolation result. The bottom row enlarges the central 
region {x,y G [—0.5, +0.5] and z G [0, 1]) outlined by the small cube in the top row. 
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Fig. 6. — CASE LLl: The history of the relaxation to force-free equihbrium. (a) Evolution 
of residual res(B) with time (or the iteration steps); (b) evolution of the maximum and 
average velocity; and (c), (d) evolution of the metrics for the central region (marked by (c)) 
and the entire volume (marked by (e)). 




Fig. 7.— Same as Figure |4] but for CASE LL2. 
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Fig. 10. — Chromospheric test case of the van Ballegooijen reference model: magnetic field 
lines with contour of on the bottom surface, (a) the reference model, (b) the extrapolation 
result. The bottom row enlarges the central region outlined by the small cube in the top 
row. 
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Fig. 11. — Chromospheric test case of the van Ballegooijen reference model: same as Fig- 
ure [T0| but projected onto the x-y planes. 
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Fig. 12. — The van Ballegooijen reference model: different views of the fluxrope. (a) 3D 
view, (b) projection on the x-y plane and (c) on the x-z plane with the white line denoting 
the bottom. Since the middle of the fluxrope lies very close to the bottom, the z scale in 
panel (c) is enlarged twice for a better view of the field lines. 
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X 



Fig. 13. — Photospheric test case of the van Ballegooijen reference model: magnetic field 
lines with contour of on the bottom surface, (a) 3D view and (b) x-y plane projection of 
the extrapolation result. The bottom row enlarges the central region outlined by the small 
cube in the top row. 
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